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Abstract 

Azimuthal asymmetries in signals of non vertical showers have been observed in 
ground arrays of water Cherenkov detectors, like Haverah Park and the Pierre Auger 
Observatory. The asymmetry in time distributions of arriving particles offers a new 
possibility for the determination of the mass composition. The dependence of this 
asymmetry on atmospheric depth shows a clear maximum at a position that is 
correlated with the primary species. In this work a novel method to determine mass 
composition based on these features of the ground signals is presented and a Monte 
Carlo study of its sensitivity is carried out. 
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1 Introduction 



The determination of the nature of ultra high energy cosmic rays (UHECR) 
is a crucial point to help understanding their origin, acceleration mechanisms 
and propagation from the sources to the Earth. At energies below 10^^ eV, 
both charge and mass can be measured directly using space detectors, how- 
ever, the properties of cosmic rays of the highest energies have to be inferred 
from the features of the shower induced in the atmosphere. Air shower experi- 
ments are either ground arrays of detectors that trigger in coincidence when a 
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shower passes through them, or optical detectors that observe the longitudinal 
development of the extensive air shower (EAS) p!|2l[3llHl5]|6] . 

The measurement of the primary mass in EAS experiments is known to be 
very difficult due to the large fluctuations resulting from the statistical nature 
of the shower development, in particular those associated to the depth and 
the number of particles produced in the first interactions. Furthermore, the 
interpretation of data to determine mass composition has to be obtained by 
comparison with Monte Carlo predictions dependent on high energy hadronic 
models. With increasing primary energy, this task becomes more difficult as 
the gap to the energy range studied in accelerator experiments increases and 
the hadronic interaction properties have to be extrapolated over a wide range. 
One of the main sources of uncertainties in any analysis to determine mass 
composition comes from the different predictions for different hadronic inter- 
action models. 




Fig. 1. Schematic view of the shower geometry. The vertical projection of the incom- 
ing direction into the shower plane (C = 0) defines two regions, "early" < C < |) 
before the shower core impact point and the opposite "late" region. Note the differ- 
ent amount of atmosphere traversed by the particles reaching the detectors in each 
region. 

The distribution of shower maximum, X^ax, that is the atmospheric depth 
at which the number of charged particles in the EAS is maximum, is sensi- 
tive to the composition of cosmic rays. Protons produce deeper showers with 
fluctuations larger than those of heavier nuclei. Therefore, for a given primary 
energy, the < X^ax > value and its fluctuations decrease with heavier primary 
mass. This is the principle in which separation methods using < X^ax >? its 
fluctuations and the elongation rate, d(Xjnax)/d(log£') [7118] . as measured by 
fluorescence detectors, are based. 

In ground array experiments the analysis is usually performed by projecting 
the signals registered by the detectors into the shower plane (see Figure [I]) 
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and thus, neglecting the further shower evolution of the late regions. As a 
consequence, for inclined showers, the circular symmetry in the signals of 
surface detectors is broken. This results in a dependence of the signal features 
on the azimuth angle in the shower plane [9], mainly due to the different 
amount of atmosphere traversed by the shower particles jlOj . 

Evidences of azimuthal asymmetries in the signal size were first observed at 
Haverah Park fllj. Recently the Pierre Auger Observatory has found in ad- 
dition, the expected asymmetry in the particle arrival time distributions [9J. 
The observation of these asymmetries for incoming directions with zenith an- 
gle smaller than 60°, has been possible at the Pierre Auger Observatory due 
to the large size of the array and the high time resolution electronic of the 
surface detector stations [12]. The design of the observatory allows measuring 
this feature of EAS which, as demonstrated below, carries very valuable infor- 
mation on to the chemical composition of cosmic rays. First results showing 
the sensitivity to primary species at the Pierre Auger Observatory have been 
presented in [13] and [H]. 

In this work it is shown that the asymmetry in risetime, ti/2, defined as the 
time to reach from 10% to 50% of the total integrated signal in each station, 
is related to the shower stage of development [T^15lll6j . Thus, for a given 
primary energy E, the asymmetry depends on zenith angle 6 of the primary 
cosmic ray in such a way that its behavior versus sec^ is reminiscent of the 
longitudinal development of the shower. This "longitudinal development of the 
asymmetry" is strongly dependent on the nature of the primary particle. The 
method presented here is quite general and, in principle, might be applied to 
other timing parameters describing the time signal structure as well as other 
shower observables like signal size which was observed to be less sensitive to 
mass composition. 

The analysis described in this work is based on Monte Carlo simulations car- 
ried out with the code aires [17] using the hadronic interaction models 
QGSJETll(03) [H] and sibyll 2.1 [19]. The aires generated showers were 
subsequently used as input in the detector simulation code and finally recon- 
structed using for both tasks the official Offline reconstruction framework of 
the Pierre Auger Observatory [2U] . 

The paper is organised as follows. In section [2] the relationship between asym- 
metry in the time structure and shower evolution is discussed in detail. A brief 
description of the Pierre Auger Observatory, which has been selected as a case 
study of this novel technique is given in section [31 The features of the Monte 
Carlo sample used for our analysis are described in section HI The core of the 
method is presented in section [5] where several mass sensitive parameters are 
defined. The procedure for the determination of the longitudinal asymmetry 
development and the definition of the discrimination parameters is firstly pre- 
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sented in I5.1[ The energy dependence of the parameters as predicted by two 
different hadronic models is discussed in I5.2[ In section [6] the relationship of 
our discrimination parameters with the shower maximum depth is analysed. 
Finally, a statistical method to estimate the primary mass composition under 
a two-component (p-Fe) assumption is presented in section [7| with a discussion 
on systematic uncertainties. 



2 Asymmetry in the time structure as an indication of shower 
evolution 

Most of the observables sensitive to composition aim at being like a snapshot 
of the shower development and thus, they are correlated with both X^^x and 
the observation depth. The time distribution of the signals recorded by the 
surface detector in EAS experiments contains implicitly the information of 
the shower development. Therefore, it is natural to expect a dependence of 
the mean value of risetime and its azimuthal asymmetry with the atmospheric 
depth traversed by the shower and thus with the zenith angle of the incoming 
cosmic ray direction. 




Fig. 2. Schematic view of shower development when arriving at three different zenith 
angles 

The sensitivity of timing parameters to primary composition can be explained 
on the basis of the dominance in the different time regions of the signal of 
the electromagnetic and muonic component. The first portion of the signal is 
dominated by the muon component which tends to arrive earlier and over a 
period of time shorter than that of the electromagnetic particles (EM), which 
are spread out on time. 

The relation between asymmetry and shower evolution is sketched in Figure 
|2l where three different scenarios are presented for a shower with a given Amax 
value, arriving at three different zenith angles. The attenuation (early-late) of 
the EM component depends on the difference in the path travelled by par- 
ticles to reach the detector stations. In case (a), i.e. vertical shower, there is 
no difference in the paths of the EM component so there is no asymmetry. 
As the zenith angle increases (case b) the difference in the attenuation of 
the EM component due to different travelled paths increases and, as a result, 
early-late asymmetry appears. At zenith angles smaller than 30° there is an 
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additional effect affecting tlie asymmetry wlien using water Clierenkov detec- 
tor for surface arrays, due to a combination of tlie geometry of tliese stations 
and the arrival direction of individual particles [THII^ . At very large zenith 
angles (case c) the EM component has been strongly absorbed before reach- 
ing the detector. Therefore, the asymmetry decreases with 6 since the larger is 
the angle the smaller is the contribution of the EM component. Note that the 
muonic component is basically asymmetry free. Then for a given zenith angle 
the asymmetry gives information of the stage of development of the shower. 
According to the above arguments, a plot of asymmetry against sec 9 is ex- 
pected to have a maximum which is correlated with the longitudinal shower 
evolution. 

The slant depth traversed by the shower particles, t', can be expressed using a 
Taylor expansion around tg = t sec6, where t is the vertical atmospheric depth 
at ground. Since the azimuth angle correction is small compared to the total 
slant depth, in a first approximation, one can keep only the first term which 
is equivalent to using a linear function in cos^ to describe asymmetries |10j. 
Thus, the atmospheric slant depth around the shower axis, can be expressed 
at fixed distance from the core, r = tq, by 

t'{C,r = To) = t sece(l + 5cosC) = ts + At,(C) (1) 

A generic time distribution for vertical showers r(r, t) depending on atmo- 
spheric depth and core distance r, becomes for inclined showers 

r{r,t) ^ T{r,t\C,e)) (2) 

where r is measured in the shower plane. 
A Taylor expansion of r around tg gives 

r(r,C) = r(r,t'(C)) = r(r,t, + At,(C)) (3) 

T{rX) = T{r,Q + ^\t^Atg{0 + ... (4) 
dr 

r(r,C) = r(r,t,) + — t^B cosC + ... (5) 

(91nr 

r(r, C) = r(r,t,)(l + ^|t, 5 cosC + ...) (6) 



Keeping only the first term of the expansion, equation ([6]) can be expressed 



as 



r(r, C) = a + 6cosC (7) 
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with 

a = T{r, t sec 9) (8) 
b dlnr , , , 

The asymmetry factor b/ a which depends on tg, can be used as a measure of the 
logarithmic variation of r with slant depth. This parameter is an indicator of 
the shower evolution and hence it provides a measure of the composition of the 
primary particle. The dependence of the asymmetry factor ^ with sec^ allows 
one to find new observables useful for determining the mass composition, as 
it will be shown in Section O 



3 The Pierre Auger Observatory 



The technique shown in this work relies in the fact that ground detectors 
located symmetrically around the shower axis are spread out over large areas. 
Therefore, the EM component of the late signal is attenuated as compared with 
that of the early one. This attenuation is only significant if the shower front is 
very extensive. The smaller is the primary energy, the shorter is the difference 
in atmospheric depth between early and late stations and thus, the weaker is 
the observed asymmetry. For instance, a reduction of the energy shower from 
10^^ to 10^"^ eV lowers the maximum distance between fired detectors by about 
a 50% and so the same percentage in the attenuation difference. On the other 
hand, the asymmetry in the time structure of the signal can be observed if 
the detectors are sensitive to both the EM and muonic components. Notice 
that the attenuation of the EM component translates into a decrease of the 
signal risetime because the earlier muon contribution is basically the same in 
the early and late regions (see Figure 2). Thus, this method can be applied 
in air-shower arrays designed for the detection of UHECR which are able to 
register information from both components. The recently completed southern 
site of the Pierre Auger Observatory [B] fulfils these requirements and therefore 
the technique is very suitable for studies of mass composition using its ground 
array. The northern site planned to be built in Colorado USA, will also be 
able to use this technique very efficiently with a larger array. The two sites 
will ensure full sky coverage. Both have been designed to use hybrid detection 
(a ground array of Cherenkov detectors overlooked by fluorescence telescopes) 
to record EASs produced by cosmic rays of energies greater than 10^^ eV. 
The Surface Detector (SD) of the southern site, situated near the town of 
Malargiie, in western Argentina, consists of 1600 stations equally spaced on 
a triangular grid (1.5 km) over an area of approximately 3000 km^. Each SD 
station is a water Cherenkov tank, 1.2 m high and top area of 10 m^ with an 
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internal reflecting coating. Three 9 inches PMTs overlook the water, and their 
signals are recorded by local digitisation electronics with a 40 MHz sampling 
rate. The southern site Fluorescence Detector (FD) consists of 4 eyes with 6 
telescopes each located at the border of the SD array overlooking it. 

The SD records the shower front, by sampling the particle density at ground 
level, with a duty cycle of 100%. The FD measures the fluorescence light 
emitted as the shower develops through the atmosphere. As it can only operate 
on clear, moonless nights, its duty cycle is about 10%. This unique combination 
of both techniques in a hybrid detector offers huge advantages in particular, 
for the determining energy [22]. With respect to primary mass estimation, 
hybrid events provide a direct measurement of X^ax- However, the bulk of 
events collected by the Observatory have information only from the surface 
array, making SD observables as the one presented in this paper very valuable 
for composition analysis at the highest energies. 



4 Monte Carlo simulated data 

The sensitivity of the method proposed here for mass composition measure- 
ments has been studied, as mentioned above, with simulated showers generated 
with AIRES 2.8.3 [H] using both hadronic interaction models QGSJETll(03) 
and SIBYLL 2.1. 

The generated data sample contains a total of 2x10^ showers, initiated by 
proton and iron nuclei. The energy and zenith angle values follows: 

• logio(^/eV) = 18.5, 19.0, 19.25, 19.5, 19.75, 20.0 

• 9 (deg) = 32, 36, 41, 45, 49, 53, 57, 60, 63 

The AIRES generated showers were subsequently used as input in the detector 
simulation code and flnally reconstructed using for both tasks the official Of- 
fline reconstruction framework of the Pierre Auger Observatory |20j. Detector 
simulation includes the generation of ground signals in the water Cherenkov 
stations of the Pierre Auger Observatory. The surface detector simulation has 
been tested and proved to be in good agreement with experimental data [23j. 

The reconstructed values of the primary energy, arrival direction and core 
location of the shower have been used in this analysis. The analysis is limited 
to the range of reconstructed 6 value above 30° and below zenith angle of 63° 
where the asymmetry effect due to the shower evolution dominates. In addition 
standard flducial cuts on SD stations have been applied including minimum 
and maximum distances to the core, signal sizes and good reconstruction of 
global shower parameters. 
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5 Mass sensitive parameters 



In this section several parameters related to the longitudinal asymmetry de- 
velopment will be defined and their sensitivity to the primary mass analysed. 

Due to the intrinsic fiuctuations in extensive air showers and the limited sam- 
pling of the shower front recorded by the SD, it is not possible to obtain the 
mass composition in a shower by shower basis. That is, instead of measuring 
the asymmetry in individual showers, the approach used in this work consists 
of using the mean value of the asymmetry factor, as defined above (equation 
[7]), for all showers in a certain interval of energy and zenith angle. 

5.1 The longitudinal asymmetry development 

The procedure used to get the average longitudinal development of the asym- 
metry, can be summarised in the following steps: 

• Select events in bins of reconstructed energy and sec^ values. For these 
events the risetime of those stations passing the corresponding cuts is de- 
termined. For each interval of E, sec 6 and (, the mean value and the stan- 
dard deviation of the corresponding ti/2/r distribution, is calculated. It is 
worth mentioning that risetime grows nearly linear with the core distance 
(see |25] and references therein) and thus the ti/2/r value is more sensitive 
for the asymmetry as the whole range of r can be used for the analysis. 

• For each [E, 9) bin, a fit of < txj^jr > to a linear cosine function of C 
(equation [7j) provides the asymmetry factor hja. Figure [3] shows as an ex- 
ample, the dependence of < ^1/2/?" > with azimuthal angle for Monte Carlo 
samples for both primaries, proton (left) and iron (right), at E = 10^^ eV 
and four zenith angles, 32°, 45°, 53° and 60°. The line represents the result 
of the fit to equation [71 From these plots it can be seen that the mean 
value of < ti/2lr > decreases with the zenith angle as expected. Besides, 
the asymmetry increases with the zenith angle up to a maximum value and 
then decreases for larger angles. 

• For each primary type and energy interval the dependence of the asymmetry 
factor on sec6' is studied. In all cases the plot b/a versus In(sec^), i.e. the 
asymmetry longitudinal development, is quite symmetric and shows a clear 
maximum. The position of the maximum is obtained by fitting a normal 
function to the values of 6/a in bins of ln(sec6') . FigurelHshows as an example 
the results for both primaries at several energies. 

This asymmetry longitudinal development can be described by means of three 
parameters, as indicated in Figure [5] (top-left): XAsymMax, the position of the 
maximum asymmetry, i.e. the sec^ value for which b/a maximises; AsymHeight, 
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Fig. 3. Risetime versus azimuth angle in the shower plane for primary energy of 
10^^ eV and different zenith angles. Proton (left) and iron (right). 






0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



Fig. 4. Asymmetry longitudinal development at primary energies 10^^'^, 10^^, 
]^0i9.25^ ]^0i9-^ 10^^-75 and 10^0 eV (top to bottom from left to right ). Monte Carlo 
proton (dashed line) and iron (solid line) are shown. 
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Fig. 5. Top-left: parameters describing the asymmetry longitudinal development. 
The other figures show the parameter dependence with primary energy for 
QGSJETii(03) and sibyll 2.1. The lines in Figure (bottom-right panel), are just to 
guide the eye. 

the height at maximum, i.e. the maximum b/a value, and XAsymWidth, the 
half width at half maximum of the Gaussian function. 



5.2 Energy dependence and hadronic models 



In Figure [5] (top-right panel and bottom panels), the values of the param- 
eters obtained from the above mentioned fit have been represented against 
primary energy. Results for both hadronic models are included in the plots. 
The error bars come from the fitting uncertainties. The bottom-left plot shows 
that XAsymMax grows linearly with \ogE. The corresponding linear fits (con- 
tinuous lines) of both primary types are clearly separated, thus allowing dis- 
crimination of heavy and light primaries. On the other hand, the AsymHeight 
parameter shows a dependence with logE which nearly follows a parabolic 
function. Results of the corresponding fit are shown for iron (red solid line) 
and proton (blue solid line) in Figure at the top-right panel. The values of 
AsymHeight for heavy and light primaries for primary energies above 10^®'^ 
eV are well separated. Finally, XAsymWidth is nearly independent on pri- 
mary energy with a value very similar for both primaries and therefore this 
parameter does not allow separation between primaries. The results seem to 
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indicate that the correlation of the asymmetry parameters with X^ax has a 
shght dependence on the hadronic interaction model. This will be addressed 
in the next section. As can be seen in Figure 0, predictions for iron do not 
show a remarkable dependence on the model. On the contrary, there is a clear 
difference for proton primaries, particularly at high energies. This is also a 
feature of the Xmax versus E plots (elongation rate) where a similar behaviour 
at high energies is observed for both models |26j . 



6 Asymmetry parameters and depth of shower msLximum 

It was mentioned before that X^ax is the main observable related to composi- 
tion in fluorescence measurements. Thus, it is desirable to study the correlation 
between our asymmetry mass sensitive parameters measured with a surface 
detector and the position of shower maximum. To this end, the steps described 
above were repeated but instead of grouping separately p and Fe events by 
primary energy, they were grouped in bins of X^ax mixing both primaries. 
The Xinax values used in these plots, are those of the simulated EAS. 

In Figure [6] the correlation with Xmax is shown for the three parameters: 
XAsymMax, AsymHeight and XAsymWidth using both hadronic interaction 
models QGSJETll(03) and sibyll 2.1. 

As expected from Figure [5], XAsymMax and AsymHeight present a strong cor- 
relation with Xmax which is nearly independent on the hadronic model. This is 
an encouraging result reaffirming that the observed azimuthal asymmetry is a 
reliable mass estimator, providing accurate models to describe EAS. Certainly, 
the correlation of the SD parameters with the position of shower maximum 
might be also useful to provide a measurement of < X^ax > using only the 
data of surface detector. 



7 Extracting the mean primary mass 

As already shown, the asymmetry longitudinal development, described by the 
three parameters above defined, is sensitive to the mass composition. The de- 
pendence of these parameters on composition has been determined as follows. 
Firstly, a two component composition (proton-iron) has been assumed. For 
each interval in E, sec 6 and (, the risetime for proton and iron showers has 
been obtained following the procedure described in section [51 The correspond- 
ing risetimes for a number of mixtures have been inferred as the average value 
for both components weighted with the corresponding composition factors. 
The percentage of iron and proton events in the sample was changed in steps 
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of 10%. As an example, the asymmetry longitudinal development curves at 
10^^ eV showers for different compositions are shown in Figure [7] (top-left). 
A smooth transition between pure proton to pure iron compositions can be 
observed. This procedure allows us to determine the dependence of the three 
discrimination parameters with composition, e.g. the Fe fraction, named in 
the following xpe- The result for 10^^ eV is shown in Figure [71 These plots 
show that XAsymMax and AsymH eight are strongly correlated with xpe while 
the correlation is much weaker for XAsymWidth, as expected from the results 
of section [51 For all energy intervals the three parameters follow a linear be- 
haviour, decreasing with the iron fraction. The grey area in Figure [71 represents 
the statistical uncertainties, determined by the errors bars from the previous 
fits (see Figures [31 and [5]) . The result for a linear fit of these three separation 
parameters is shown in this figure (solid line). The parameterisation of this 
dependence will be used later to extract the primary mass. The composition 
of a sample, i.e. the xpe value (either a Monte Carlo test sample or real data) 
can be determined by minimising the following function 
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where a, h and m stand for the parameters XAsymWidth, AsymHeight and 
XAsymMax, respectively; the subindex s indicates the corresponding values of 
the sample to be studied. Aa, Ah and Am represent the standard deviation of 
the separation parameters. The functions fi{xpe) stand for the parameterisa- 
tion of the dependence of each parameter with the iron fraction xpe obtained 
in Figure [71 The grey area indicates the uncertainty in the parameterisation 
due to limited Monte Carlo statistics, and it will be taken into account as a 
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Fig. 7. Asymmetry development for the different samples with mixed composition, 
going from pure proton to pure iron in steps of 10%. Position of the maximum, 
sigma and height of the longitudinal profile as a function of the mixture fraction. 
1 corresponds to 100 % iron and to 100 % proton. Parameters shown here as an 
example correspond to primary energy of 10^^ eV. 



systematic error. 

To check the rehability of this technique we have reconstructed the xpe value 
of several samples containing a known fraction of proton and iron showers 
from our Monte Carlo data by minimising A^. This test has been carried 
out using mixtures of showers containing the same events used previously to 
get the results on Figure [71 The results are shown in Figure [HI for three 
different samples: 100% iron, i.e. Xpe = 1, (upper curve), 50% iron and 50% 
proton, i.e. xpe = 0.5, (middle curve) and 100% proton, i.e. xpe = 0, (lower 
curve). Deviations of reconstructed xpe values from the real ones are below 
0.04 (i.e. errors in the Fe percentage below 4 units). The error bars represent 
the uncertainty due to our limited statistics and amount about Axpe = 0.1. 
This of course can be reduced using larger MC samples. 

A study of systematic errors has been performed. Firstly, the full analysis has 
been carried out using both hadronic models, using the parameterisations from 
QGSJETll(03) and then applying that to sibyll 2.1 data set. The prediction 
in the xpe value for both models are different by an amount which ranges up 
to Axpe= 0.14. 
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Fig. 8. Consistency check, iron fraction with three different samples: 100% 
iron(upper curve), 50% iron and 50% proton (middle curve), 100% proton(lower 
curve) 



Due to the influence of the models the primary energy reconstructed for SD 
events using Monte Carlo, deviates from hybrid reconstructed primary ener- 
gies. The effect of this deviation in our method has been studied by comparing 
the composition obtained using both the MC E value (i.e the input for the 
Monte Carlo generation) and the one reconstructed with the Auger analysis 
software. Differences in xpe lower than 0.04 have been found. 



A possible contribution due to the cuts applied in our method (see section HI) 
has been studied by comparing results with different cuts. The corresponding 
contribution results in Axpe = 0.03. 



It should be mentioned at this point that taking the estimated event rate 
of the Pierre Auger Observatory, a similar resolution than the one presented 
in Figure [8] can be achieved at present for energies below 10^^ eV. With 
increasing statistics collected by the full array in the near future, the shown 
resolution will be achieved at higher energies. 
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8 Conclusions 



A novel method to determine the mass composition of primary cosmic rays has 
been developed using the azimuthal asymmetry in arrival time distribution of 
secondary particles at a given observing level. The approach relies on statistical 
grounds and thus provides a mean mass composition of a set of showers at a 
given energy. 

The main idea behind the method is to reconstruct a longitudinal develop- 
ment of the observed asymmetry which is reminiscent of the longitudinal de- 
velopment of the extensive air shower. A detailed analysis using the risetime 
of the signal in water Cherenkov detectors for the case of the Pierre Auger 
Observatory was presented. It was shown that both the atmospheric depth 
corresponding to the position of maximum asymmetry and the value of the 
maximum asymmetry are sensitive to primary mass. These parameters mea- 
sured by the surface detectors were shown to correlate with the position of 
shower maximum, X^ax- 

The method was validated using hypothetical data samples corresponding to 
pure proton, pure iron and a mixed composition. Systematic uncertainties 
affecting the determination of primary composition were investigated. As ex- 
pected, the dominant source of uncertainties comes from the lack of knowledge 
of hadronic interaction models, which amounts to < 14% out of a total of 18% 
in the estimation of Fe fraction. The analysis indicates that at the event rate 
collected by the Pierre Auger Observatory, a very good separation of heavy 
and light elements can be achieved with present data, for energies below 10^^ 
eV, while at higher energies at least one order of magnitude more data would 
be needed to achieve the shown resolution in this paper. 
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